#### Lab 8 #####

### Set up
daphnia = read.table("V:\\My Documents\\ENT420\\lab\\lab6_data.txt", sep="\t", header=T)

### Problem 1

## a.  rename this to p
p = length(daphnia[daphnia$past==1,]$past) / length(daphnia$past)

##b. write it out

##c.
therebefore = daphnia[daphnia$past==1,]
lost = therebefore[therebefore$now==0,]
prob = 1 - exp(-0.5)
dbinom(length(lost$now), length(therebefore$past), prob)
# 0.1903814

##d.
probdist = dbinom(length(lost$now), length(therebefore$past),seq(0,1,by=0.01))
maxlikelihood = max(probdist)
seq(0,1,by=0.01)[probdist==max(probdist)] # .31, extinction probability
E = -log(1-.31)
#plot(probdist)

##e.
unocc = daphnia[daphnia$past==0,]
numcol = unocc[unocc$now==1,]
mean(unocc$now) 	# 0.18 - proportion now occupied
length(unocc$x) 	# 72
length(numcol$x)	# 13


##f.
#
dbinom(length(numcol$x),length(unocc$x),(1-exp(-0.5*p)))
prolike = dbinom(length(numcol$x),length(unocc$x),seq(0,1,by=0.01))
#likelyC = max(prolike) # 0.1214
likelyC = seq(0,1,by=0.01)[prolike==max(prolike)] 
C = -log(1-likelyC) / p
C 
#plot(prolike)


##h.
(C - E) / C 		# predicted equilibrium occupancy

#glm(now ~ past, family = binomial(), data = daphnia)



##### Problem 2 #####
#dim() gives the numbers of locations (rows), p, and numbers of variables (cols) in a
#data.frame. The distance matrix will need to be a p x p matrix. Let's therefore pick the
#number of locations and initiate the p x p distance matrix (of missing values):
p=dim(daphnia)[1] #nos of rows in the daphnia data-set
dmat = matrix(NA,nrow=p, ncol=p) #initiating the p x p matrix
#calculating the distances and filling them into the correct place in the dmat (remember:
#indents are ignored by R, but help in our coding)
for(i in 1:p){
for(j in 1:p){
dmat[i, j] = sqrt((daphnia$x[i]-daphnia$x[j])^2+(daphnia$y[i]-daphnia$y[j])^2)
}
}
#Advanced - double loops are computationally expensive so R has implemented a
#shortcut
#dmat = sqrt(outer(daphnia$x, daphnia$x, "-")^2+outer(daphnia$y,daphnia$y,"-")^2)
#To study how presence/absence depends on connectivity (and area) we may FIRST
#postulate ('guess') a value for a at around 100km. (so as to get our code together).
# if a is 100
a=100
#a=33.24
#exponentially weighted matrix (given a = 100) weighed by past occupancy
wmat = exp(-dmat/a)*daphnia$past
#(assumed) proportional immigration of target patches (given a) the summed
#contribution across the donors (non-donors having been nulled out) is:
connectivity=apply(wmat, 2, sum)
#logistic regression of presrence as a function of area and connectivity (given that a is
#100). Obviously, previously occupied patches are extremely connected since local
#persistence can be secured through non-extinction of that local lake.
Fit=glm(now~log(area)+ connectivity, data=daphnia, family=binomial())
#the associated log-likelihood is -1/2*deviance
-1/2*Fit$deviance  # for a==100, -39.62

a=seq(1, 400, length=100) #candidate a's
llik = rep(NA, 100) #initiating a vector of NA's to hold
#the corresponding log-likelihoods
for(i in 1:100){ #loop across the a's
#calculate corresponding proportional immigration summed
#across donor pops (non-donors having been nulled out)
connectivity =apply(exp(-dmat/a[i])*daphnia$past, 2, sum)
fit=glm(now~log(area)+ connectivity, data=daphnia,
family=binomial())
#the log-likelihood is -1/2*deviance, to be stored in the llik vector
llik[i] = -1/2*fit$deviance
}
#to look for evidence of connectivity effects we inspect how the log-likelihood
#is a function candidate a's

plot(a, llik, log="x", type="b")


#a.)
max(llik) # -38.46
a[llik==max(llik)]  #33.24

#b.)
a=33.24
wmat = exp(-dmat/a)*daphnia$past
connectivity=apply(wmat, 2, sum)
Fit=glm(now~log(area)+ connectivity, data=daphnia, family=binomial())
-1/2*Fit$deviance  # -38.46018

#c.) 
symbols(daphnia$x,daphnia$y,circles=log(daphnia$area), inches=.2)